#importing libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
from tqdm import tqdm_notebook
from igraph import Graph
import igraph
import Neuroscience_Project_modules as npm
import bct
import random
##paths
open_eyes = "./data/S004R01.edf"
closed_eyes = "./data/S004R02.edf"
freq_alpha = 10 ###10Hz Alpha_Wave_Frequency
Estimate functional brain connectivity among 64 channels using one of the MVAR estimators: Partial Directed Coherence (PDC), Direct Transfer Function (DTF). Select one relevant frequency value. Apply a threshold so that the resulting binary connectivity matrices have network density equal to 20%. Create a graphical representation of the binary adjacency matrix.
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(PDC_open.binary_adjacency_matrix)
plt.imshow(PDC_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(PDC_closed.binary_adjacency_matrix)
plt.imshow(PDC_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Perform task 1.1 using both estimators (PDC and DTF).
DTF_open = npm.Connectivity_Graph('DTF')
DTF_open.import_data(open_eyes)
DTF_open.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(DTF_open.binary_adjacency_matrix)
plt.imshow(DTF_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
DTF_closed = npm.Connectivity_Graph('DTF')
DTF_closed.import_data(closed_eyes)
DTF_closed.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
print(DTF_closed.binary_adjacency_matrix)
plt.imshow(DTF_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
#Models Order computed from Akaike criterion
order_open = 7
order_closed = 6
Perform task 1.1 using thresholds yielding the following density values: 1%, 5%, 10%, 20%, 30%, 50%.
m_open = []
m_closed = []
densities = [0.01,0.05,0.1,0.2,0.3,0.5]
for d in tqdm_notebook(densities):
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=d,order=order_open)
m_open.append(PDC_open.binary_adjacency_matrix)
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=d,order=order_closed)
m_closed.append(PDC_closed.binary_adjacency_matrix)
m = [m_open,m_closed]
cols = ['Case: {}'.format(col) for col in ['open','closed']]
rows = ['Density: {}'.format(row) for row in densities]
fig, axes = plt.subplots(nrows=6, ncols=2, figsize=(20,12))
for ax, col in zip(axes[0], cols):
ax.set_title(col)
for ax, row in zip(axes[:,0], rows):
ax.set_ylabel(row, rotation=90, size='large')
r = 0
for row in axes:
c = 0
for col in row:
col.imshow(m[c][r],cmap='Greys',interpolation='nearest')
c += 1
r += 1
fig.tight_layout()
fig.subplots_adjust(left=0.05, top=0.95)
plt.show()
Considering the subset of 19 channels suggested in Figure 1 and Table 2, estimate the connectivity using PDC or DTF and apply a statistical validation method (asymptotic statistics 7 , resampling procedure 8 ,...) to filter out values that are not significantly different from 0 (PDC(i, j) ≠ 0 with p < 5%).
channels = 'Fp1 Fp2 F7 F3 Fz F4 F8 T7 C3 Cz C4 T8 P7 P3 Pz P4 P8 O1 O2'.split(' ')
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes,channels=channels)
PDC_open.connectivity(freq=freq_alpha,order=order_open)
PDC_open.significance(signf_threshold=0.05,channels=channels,visual=False,path=open_eyes,Nrep=1000)
#PDC_open.significance_matrix
plt.imshow(PDC_open.significance_matrix,cmap='jet_r')
## associated binary adjacency matrix after insignificant connections removal
plt.imshow(PDC_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes,channels=channels)
PDC_closed.connectivity(freq=freq_alpha,order=order_closed)
PDC_closed.significance(signf_threshold=0.05,channels=channels,visual=False,path=open_eyes,Nrep=1000)
#PDC_closed.significance_matrix
plt.imshow(PDC_closed.significance_matrix,cmap='jet_r')
## associated binary adjacency matrix after insignificant connections removal
plt.imshow(PDC_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Make a topographical representation of the networks (see example in Figure 2). Cartesian coordinates of planar representation of EEG channels are available in Table 3 (see also the file channel_locations.txt).
(the choice of this task is advised in the case of 19-channel networks and/or density ≤ 5%).
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=0.05,order=order_open)
PDC_open.show_graph('PDC_case_open')
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=0.05,order=order_closed)
PDC_closed.show_graph('PDC_case_closed')
Perform task 1.1 considering a second frequency value belonging to a different EEG rhythm with respect to the first choice.
freq_beta = 20 ###20Hz Beta_Wave_Frequency
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_beta,threshold=0.2,order=order_open)
print(PDC_open.binary_adjacency_matrix)
plt.imshow(PDC_open.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_beta,threshold=0.2,order=order_closed)
print(PDC_closed.binary_adjacency_matrix)
plt.imshow(PDC_closed.binary_adjacency_matrix,cmap='Greys',interpolation='nearest')
Compute binary global (clustering coefficient, path length) and local (degree, in/out-degree) graph indices. List the highest 10 channels for local indices.
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Search in the literature a definition of small-worldness index (i.e. an index describing the small-world organization of a network) and compute it.
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
C = IDX.global_clustering_coefficient(False)
L = IDX.average_path_length(False)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nGraph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L '+'\x1b[0m' + '= ' + str(L),'\x1b[1;39m'+ ' C '+'\x1b[0m' + '= ' + str(C) )
APL_rnd_list = []
GCC_rnd_list = []
for cnt in range(10):
G_r_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_r_Adj = bct.randmio_dir(G_r_Adj, 10)[0]
G_r_Adj = [list(r) for r in G_r_Adj]
G_r_Adj = [[G_r_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_r = Graph.Adjacency(G_r_Adj)
APL_rnd_list.append(npm.average_path_length(G_r, False))
GCC_rnd_list.append(npm.global_clustering_coefficient(G_r, False))
L_r = np.mean(APL_rnd_list)
C_r = np.mean(GCC_rnd_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Random Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_r '+'\x1b[0m' + '= ' + str(L_r),'\x1b[1;39m'+ ' C_r '+'\x1b[0m' + '= ' + str(C_r) )
APL_ltc_list = []
GCC_ltc_list = []
for cnt in range(10):
G_l_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_l_Adj = bct.latmio_dir(G_l_Adj, 10)[0]
G_l_Adj = [list(r) for r in G_l_Adj]
G_l_Adj = [[G_l_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_l = Graph.Adjacency(G_l_Adj)
APL_ltc_list.append(npm.average_path_length(G_l, False))
GCC_ltc_list.append(npm.global_clustering_coefficient(G_l, False))
L_l = np.mean(APL_ltc_list)
C_l = np.mean(GCC_ltc_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Lattice Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_l '+'\x1b[0m' + '= ' + str(L_l),'\x1b[1;39m'+ ' C_l '+'\x1b[0m' + '= ' + str(C_l) )
SWI = (L - L_l)/(L_r - L_l) * (C - C_r)/(C_l - C_r)
print('\x1b[1;35m'+'\nSmall World Index\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'SWI '+'\x1b[0m' + '=', SWI)
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
C = IDX.global_clustering_coefficient(False)
L = IDX.average_path_length(False)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nGraph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L '+'\x1b[0m' + '= ' + str(L),'\x1b[1;39m'+ ' C '+'\x1b[0m' + '= ' + str(C) )
APL_rnd_list = []
GCC_rnd_list = []
for cnt in range(10):
G_r_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_r_Adj = bct.randmio_dir(G_r_Adj, 10)[0]
G_r_Adj = [list(r) for r in G_r_Adj]
G_r_Adj = [[G_r_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_r = Graph.Adjacency(G_r_Adj)
APL_rnd_list.append(npm.average_path_length(G_r, False))
GCC_rnd_list.append(npm.global_clustering_coefficient(G_r, False))
L_r = np.mean(APL_rnd_list)
C_r = np.mean(GCC_rnd_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Random Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_r '+'\x1b[0m' + '= ' + str(L_r),'\x1b[1;39m'+ ' C_r '+'\x1b[0m' + '= ' + str(C_r) )
APL_ltc_list = []
GCC_ltc_list = []
for cnt in range(10):
G_l_Adj = np.array([r for r in IDX.G.get_adjacency()])
G_l_Adj = bct.latmio_dir(G_l_Adj, 10)[0]
G_l_Adj = [list(r) for r in G_l_Adj]
G_l_Adj = [[G_l_Adj[i][j].item() for j in range(len(G_r_Adj))] for i in range(len(G_r_Adj))]
G_l = Graph.Adjacency(G_l_Adj)
APL_ltc_list.append(npm.average_path_length(G_l, False))
GCC_ltc_list.append(npm.global_clustering_coefficient(G_l, False))
L_l = np.mean(APL_ltc_list)
C_l = np.mean(GCC_ltc_list)
print('\x1b[1;35m'+'\nAverage Path Length and Global Clustering Coefficient for\nEquivalent Lattice Graph\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'L_l '+'\x1b[0m' + '= ' + str(L_l),'\x1b[1;39m'+ ' C_l '+'\x1b[0m' + '= ' + str(C_l) )
SWI = (L - L_l)/(L_r - L_l) * (C - C_r)/(C_l - C_r)
print('\x1b[1;35m'+'\nSmall World Index\n' + '\x1b[0m')
print('\x1b[1;39m'+ 'SWI '+'\x1b[0m' + '=', SWI)
Compare the global indices extracted from PDC and DTF connectivity estimations.
IDX = npm.Graph_Theory_Indices("DTF")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
IDX = npm.Graph_Theory_Indices("DTF")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Study the behaviour of global graph indices in function of network density (see point 2.3 for density values).
(the choice of this task is advised in the case of selection of task 1.3).
random.seed (7890)
APL_list = []
GCC_list = []
densities = [0.01,0.05,0.1,0.2,0.3,0.5]
for d in tqdm_notebook(densities):
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=d, plot=False,order=order_open)
APL_list.append(IDX.average_path_length(False))
GCC_list.append(IDX.global_clustering_coefficient(False))
plt.plot(densities, APL_list)
plt.plot(densities, GCC_list)
APL_list = []
GCC_list = []
densities = [0.01,0.05,0.1,0.2,0.3,0.5]
for d in tqdm_notebook(densities):
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=d, plot=False,order=order_closed)
APL_list.append(IDX.average_path_length(False))
GCC_list.append(IDX.global_clustering_coefficient(False))
plt.plot(densities, APL_list)
plt.plot(densities, GCC_list)
Make a topographical representation of local indices.
def plot_graph(G, name, degree_list):
G.es['arrow_size'] = [0.1 for edge in G.es]
visual_style = {}
visual_style["vertex_size"] = degree_list*10
visual_style["vertex_color"] = "red"
visual_style["vertex_label"] = G.vs["label"]
visual_style["edge_width"] = 0.2
visual_style["layout"] = G.vs["coords"]
graph = igraph.plot(G, bbox=(0, 0, 600, 600), **visual_style)
graph.save(name + '.png')
return(graph)
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
degree_list = list(IDX.local_indegree().values())
plot_graph(IDX.G, "EO_indegree_graph", degree_list)
degree_list = list(IDX.local_outdegree().values())
plot_graph(IDX.G, "EO_outdegree_graph", degree_list)
degree_list = list(IDX.local_degree().values())
plot_graph(IDX.G, "EO_degree_graph", degree_list)
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_closed)
degree_list = list(IDX.local_indegree().values())
plot_graph(IDX.G, "EC_indegree_graph", degree_list)
degree_list = list(IDX.local_outdegree().values())
plot_graph(IDX.G, "EC_outdegree_graph", degree_list)
degree_list = list(IDX.local_degree().values())
plot_graph(IDX.G, "EC_degree_graph", degree_list)
Compare the networks obtained with the analysis 1.6 in terms of graph indices.
(the choice of this task is advised only in the case of selection of task 1.6)
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_beta, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_beta, threshold=0.2, plot=False,order=order_closed)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/o weights): "+'\x1b[0m', IDX.average_path_length(False))
print('\x1b[1;39m'+"GCC (w/o weights): "+'\x1b[0m', IDX.global_clustering_coefficient(False))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree())
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_outdegree())
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_degree())
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree()
local_outdegree_dict = IDX.local_outdegree()
local_degree_dict = IDX.local_degree()
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Perform point 2.1 considering the weighted version of the graph indices definitions.
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(open_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/ weights): "+'\x1b[0m', IDX.average_path_length(True))
print('\x1b[1;39m'+"GCC (w/ weights): "+'\x1b[0m', IDX.global_clustering_coefficient(True))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree(True)
local_outdegree_dict = IDX.local_outdegree(True)
local_degree_dict = IDX.local_degree(True)
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
IDX = npm.Graph_Theory_Indices("PDC")
IDX.import_data(closed_eyes)
IDX.connectivity(freq=freq_alpha, threshold=0.2, plot=False,order=order_open)
print('\x1b[1;35m'+'\nGlobal indices'+'\x1b[0m \n')
print("APL (Average Path Length) and GCC (Global Clustering Coefficient)\n")
print('\x1b[1;39m'+"APL (w/ weights): "+'\x1b[0m', IDX.average_path_length(True))
print('\x1b[1;39m'+"GCC (w/ weights): "+'\x1b[0m', IDX.global_clustering_coefficient(True))
print("\n")
print('\x1b[1;35m'+'Local indices'+'\x1b[0m \n')
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local outdegree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;39m'+"Local degree: "+'\x1b[0m', IDX.local_indegree(True))
print('\x1b[1;35m'+'Top 10'+'\x1b[0m \n')
local_indegree_dict = IDX.local_indegree(True)
local_outdegree_dict = IDX.local_outdegree(True)
local_degree_dict = IDX.local_degree(True)
print('\x1b[1;39m'+"Local indegree: "+'\x1b[0m\n')
print(npm.top10(local_indegree_dict))
print('\n\x1b[1;39m'+"Local outdegree: "+'\x1b[0m\n')
print(npm.top10(local_outdegree_dict))
print('\n\x1b[1;39m'+"Local degree: "+'\x1b[0m\n')
print(npm.top10(local_degree_dict))
Perform motifs analysis to investigate the presence of 3-node configurations in the networks: determine their frequency and statistical significance (motifs, anti-motifs).
## alpha waves
## threshold: 0.05
## method: PDC
## method for p: AIC
freq=10
threshold=0.05
method = 'PDC'
# Create EO connectivity graph for 10Hz frequence
# then save the edges in txt file for mfinder app
EO = npm.Graph_for_motifs(method, 'EO', density_threshold=threshold, freq=freq, folder='data/')
EO.connectivity(freq=freq,threshold=threshold,plot=True)
EO.preprocess_for_mfinder()
# performing motifs analysis with mfinder:
# mfinder ../data/EO_10Hz_0.05dt_Motifs.txt -s 3 -r 2000 -nu
# mfinder ../data/EC_10Hz_0.05dt_Motifs.txt -s 3 -r 2000 -nu
with open('data/EO_10Hz_0.05dt_Motifs_OUT.txt', 'r') as f:
print (f.read())
===============
2 motifs detected:
| ID | 6, 108 |

===============
1 anti-motifs detected:
| ID | x |

# Create EC connectivity graph for 10Hz frequence
# then save the edges in txt file for mfinder app
EC = npm.Graph_for_motifs(method, 'EC', density_threshold=threshold, freq=freq_alpha, folder='data/')
EC.connectivity(freq=freq,threshold=threshold,plot=True)
EC.preprocess_for_mfinder()
with open('data/EC_10Hz_0.05dt_Motifs_OUT.txt', 'r') as f:
print (f.read())
===============
4 motifs detected:
| ID | 38, 46, 108, 238|

===============
x anti-motifs detected

| ID | x |
For the motif with pattern A → B ← C, create a topographical representation of the networks considering only the connections involved in this configuration.
# Motif file structure: EO
with open('data/EO_10Hz_0.05dt_Motifs_MEMBERS.txt', 'r') as f:
for _ in range(9):
print(f.readline())
print('\t...')
EO.plot_motifs_36('data/EO_10Hz_0.05dt_Motifs_MEMBERS.txt',
name='data/EO_0.05_motif36', save=False)
# Motif file structure: EC
with open('data/EC_10Hz_0.05dt_Motifs_MEMBERS.txt', 'r') as f:
for _ in range(9):
print(f.readline())
print('\t...')
EC.plot_motifs_36('data/EC_10Hz_0.05dt_Motifs_MEMBERS.txt',
name='data/EC_0.05_motif36', save=False)
Choose a channel selected in parieto-occipital scalp region and determine the motifs which involve it.
df_prova = EO.find_motifs(channel='Po3')
df_prova
EC.find_motifs('Po3')
# channel selected: Po3 (57)
# in our notation (from 0 to 64) Po3 correspond to 56
# we can directly plot it with our figure
# remaking member analysis
# plot for EO (pattern A->B<-C)
EO.plot_motifs_36('data/EO_10Hz_0.05dt_Motifs_MEMBERS.txt', 'EO_motifsPlot_Po3',
single_channel=56, motifcolor="red", save=False, tipsize=0.5)
# plot for EC (pattern A->B<-C)
EC.plot_motifs_36('data/EC_10Hz_0.05dt_Motifs_MEMBERS.txt', 'EC_motifsPlot_Po3',
single_channel=56, motifcolor="red", save=False, tipsize=0.5)
Perform the same analysis described in task 3.1 considering 4-node motifs
The same analysis has been performed using mfinder. The netork used was the same of point 3.1
with open('data/EO_10Hz_0.05dt_4_Motifs_OUT.txt', 'r') as f:
for _ in range(51):
print(f.readline())
print('...\n...\nCheck EO_10Hz_0.05dt_4_Motifs_OUT.txt for more info')
with open('data/EC_10Hz_0.05dt_4_Motifs_OUT.txt', 'r') as f:
for _ in range(51):
print(f.readline())
print('...\n...\nCheck EC_10Hz_0.05dt_4_Motifs_OUT.txt for more info')
Determine number and composition (i.e. list of nodes) of the communities obtained applying one of the algorithms introduced during the course.
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=0.2,plot=True)
Make a graphical representation of the community structure in both rest conditions.
comms_mod_open,graph = PDC_open.optimal_modularity_community_detection(name='optimal_modularity_open_eyes')
graph
comms_mod_closed,graph = PDC_closed.optimal_modularity_community_detection(name='optimal_modularity_closed_eyes')
graph
Compare the community structure obtained by means of two different methods (modularity-based vs information theory-based approaches).
comms_info_open,graph = PDC_open.infomap_community_detection(name='information_theory_open_eyes')
graph
M = np.zeros([len(comms_mod_open),len(comms_info_open)])
for i in range(len(comms_mod_open)):
for j in range(len(comms_info_open)):
#JACCARD SIMILARITY
M[i,j] = len(set(comms_mod_open[i]).intersection(set(comms_info_open[j])))/len(set(comms_mod_open[i]).union(set(comms_info_open[j])))
M = pd.DataFrame(M)
plt.figure(figsize=(9,9))
ax = sns.heatmap(M,xticklabels=M.columns+1,yticklabels=M.columns+1,annot=True,cmap='viridis',linewidths=1,square=True)
ax.invert_yaxis()
plt.xlabel('Infomap',fontsize=12)
xlabels = comms_info_open.keys()
ylabels = comms_mod_open.keys()
plt.xticks(np.arange(len(comms_info_open))+0.5, xlabels,)
plt.yticks(np.arange(len(comms_mod_open))+0.5, ylabels)
plt.ylim(0,len(comms_mod_open))
plt.ylabel('Optimal_Modularity_Method',fontsize=12)
plt.title('Jaccard Similarity',fontsize=15)
plt.tight_layout()
plt.show()
comms_info_closed,graph = PDC_closed.infomap_community_detection(name='information_theory_closed_eyes')
graph
M = np.zeros([len(comms_mod_closed),len(comms_info_closed)])
for i in range(len(comms_mod_closed)):
for j in range(len(comms_info_closed)):
#JACCARD SIMILARITY
M[i,j] = len(set(comms_mod_closed[i]).intersection(set(comms_info_closed[j])))/len(set(comms_mod_closed[i]).union(set(comms_info_closed[j])))
M = pd.DataFrame(M)
plt.figure(figsize=(9,9))
ax = sns.heatmap(M,xticklabels=M.columns+1,yticklabels=M.columns+1,annot=True,cmap='viridis',linewidths=1,square=True)
ax.invert_yaxis()
plt.xlabel('Infomap',fontsize=12)
xlabels = comms_info_closed.keys()
ylabels = comms_mod_closed.keys()
plt.xticks(np.arange(len(comms_info_closed))+0.5, xlabels,)
plt.yticks(np.arange(len(comms_mod_closed))+0.5, ylabels)
plt.ylim(0,len(comms_mod_closed))
plt.ylabel('Optimal_Modularity_Method',fontsize=12)
plt.title('Jaccard Similarity',fontsize=15)
plt.tight_layout()
plt.show()
##densities threshold
densities = [0.05,0.1,0.2,0.3,0.5]
m1 = []
m2 = []
m3 = []
m4 = []
for d in tqdm_notebook(densities):
PDC_open = npm.Connectivity_Graph('PDC')
PDC_open.import_data(open_eyes)
PDC_open.connectivity(freq=freq_alpha,threshold=d,plot=False,order=order_open)
PDC_closed = npm.Connectivity_Graph('PDC')
PDC_closed.import_data(closed_eyes)
PDC_closed.connectivity(freq=freq_alpha,threshold=d,plot=False,order=order_closed)
comms_open_mod = PDC_open.optimal_modularity_community_detection(visual = False)
comms_closed_mod = PDC_closed.optimal_modularity_community_detection(visual = False)
comms_open_info = PDC_open.infomap_community_detection(visual = False)
comms_closed_info = PDC_closed.infomap_community_detection(visual = False)
m1.append(len(comms_open_mod))
m2.append(len(comms_open_info))
m3.append(len(comms_closed_mod))
m4.append(len(comms_closed_info))
# set width of bar
barWidth = 0.25
# set height of bar
bars1 = m1
bars2 = m2
# Set position of bar on X axis
r1 = np.arange(len(bars1))
r2 = [x + barWidth for x in r1]
# Make the plot
plt.bar(r1, bars1, color='#7f6d5f', width=barWidth, edgecolor='white', label='Louvain')
plt.bar(r2, bars2, color='#557f2d', width=barWidth, edgecolor='white', label='Infomap')
# Add xticks on the middle of the group bars
plt.xlabel('graph density', fontweight='bold')
plt.xticks([r + barWidth for r in range(len(bars1))], densities)
# Create legend & Show graphic
plt.title('Number of communities - Open Eyes Case')
plt.legend()
plt.show()
# set width of bar
barWidth = 0.25
# set height of bar
bars1 = m3
bars2 = m4
# Set position of bar on X axis
r1 = np.arange(len(bars1))
r2 = [x + barWidth for x in r1]
# Make the plot
plt.bar(r1, bars1, color='#7f6d5f', width=barWidth, edgecolor='white', label='Louvain')
plt.bar(r2, bars2, color='#557f2d', width=barWidth, edgecolor='white', label='Infomap')
# Add xticks on the middle of the group bars
plt.xlabel('graph density', fontweight='bold')
plt.xticks([r + barWidth for r in range(len(bars1))], densities)
# Create legend & Show graphic
plt.title('Number of communities - Closed Eyes Case')
plt.legend()
plt.show()